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Abstract 

ABS methods are a large class of methods, based upon the Egervary rank 
reducing algebraic process, first introduced in 1984 by AbafFy, Broyden and 
Spedicato for solving linear algebraic systems, and later extended to nonlin- 
ear algebraic equations, to optimization problems and other fields; software 
based upon ABS methods is now under development. Current ABS literature 
consists of about 400 papers. ABS methods provide a unification of several 
classes of classical algorithms and more efficient new solvers for a number of 
problems. In this paper we review ABS methods for linear systems and opti- 
mization, from both the point of view of theory and the numerical performance 
of ABSPACK. 

Key Words ABS methods, linear algebraic systems, feasible direction 
methods, simplex method, Diophantine equations, ABSPACK. 

AMS classification 15A03, 15A06, 65F05, 65F30 



1 Introduction 

ABS algorithms were introduced by Abaffy, Broyden and Spedicato (1984), to solve 
linear equations first in the form of the basic ABS class, later generalized as the scaled 
ABS class. They were then applied to linear least squares, nonlinear equations and 
optimization problems, see e.g. the monographs by AbafFy and Spedicato (1989) 
and Feng et al. (1998), or the bibhography by Spedicato et al. (2000) listing 350 
ABS papers. In this paper we review some of the main results obtained in the field 
of ABS methods and wc provide some results on the performance of ABSPACK, a 
FORTRAN package based on ABS methods presently under development. 

The main results obtained in almost twenty years of research in the ABS field 

can be summarized as follows: 

• ABS methods provide a unification of the field of finitely terminating methods 
for linear systems and for feasible direction methods for linearly constrained 
optimization, see sections 2 and 6; due to the several alternative formulations 
of their algebra they lead to different computational implementations of the 
algorithms, each one with a special feature that may be of advantage 

• ABS methods provide some new methods that are better than classical ones 
under several respects. For instance the implicit LX method requires the same 
number of multiplications as Gaussian elimination (which is optimal under 
very general conditions) but requires less memory and does not need pivoting 
in general; moreover its application to the simplex method has not only a lower 
memory requirement, but is cheaper in multiplications up to one order with 
respect to e.g. the Forrest- Goldfarb implementation of the simplex method via 
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the LU factorization, see section 5. Also there are ABS methods for nonhnear 
problems that are better than Newton method in terms of convergence speed 
(some ABS algorithms keep quadratic rate of convergence even if the Jacobian 
is singular at the solution) or in terms of required information (about 
Jacobian components against the full n^) 

• on some classes of significant problems ABSPACK codes have a better per- 
formance in both accuracy and speed that the codes in LAPACK (usually 
considered the best package now on the market); see section 7. 

• ABS methods have allowed to solve open problems in the literature (e.g. the 
explicit determination of Quasi-Newton updates for the sparse symmetric case, 
see Spedicato and Zhao (1993)) 

• for linear Diophantine equations ABS methods have led to a new solution 
method, that also provides the first generalization of the classical existence 
theorem of Euclides, Diophantus and Euler from a single n-dimensional equa- 
tion to a general m-equations system in n variables, see section 4. 

2 The scaled ABS class: definition and main 
properties 

Let us consider the following determined or underdetermined linear system, where 
rank{A) is arbitrary and — (ai, . . . , a^) 

Ax^b xeEJ", beR"", m<n (1) 

or 

ajx — bi = 0, i — l,...,m (2) 

The above system can be solved by the following scaled ABS class of algorithms, 
which can be shown to provide a complete realization of the general class of methods 
implementing the Petrov-Galerkin criterion (stating that the residual vector at the 
i-th iteration is orthogonal to a given arbitrary set of i — 1 linearly independent 
vectors) . 

The scaled ABS algorithm 

(A) Give Xi e arbitrary. Hi e i?"'" nonsingular arbitrary, vi e RJ^ arbitrary 
nonzero. Set i — 1. 

(B) Compute the residual = Axi — b. If = stop {xi solves the problem.) 
Otherwise compute Si — HiA^Vi. If Si ^ 0, then go to (C). If Sj = and 
T = vjri = 0, then set Xj+i = Xi, i^j+i = Hi (the i-th equation is redundant) 
and go to (F). Otherwise stop (the system has no solution). 
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(C) Compute the search vector pi by 

Pi = Hjzi (3) 
where Zi e i?" is arbitrary save for the condition 

vjAHfzi + (4) 

(D) Update the estimate of the solution by 

Xj+i = - ftiPi, ai^vjrijvj Api. (5) 

(E) Update the matrix Hi by 

//i+i ^Ei- HiA^ViwjHi/wjHiA'^Vi (6) 
where Wi G i?" is arbitrary save for the condition 

wjHiA^Vi ^ 0. (7) 

(F) If i = m, stop {xm+i solves the system). Otherwise give G R"^ arbitrary 

linearly independent from f i, . . . , f j. Increment i by one and go to (B). 

Matrices Hi, which are generalizations of (oblique) projection matrices, have been 
named Ahaffians at the First International Conference on ABS methods (Luoyang, 
China, 1991). However we must recall that these matrices were used in several little 
known papers by Egervary, predating the ABS algorithms. There are alternative 
formulations of the scaled ABS algorithms, e.g. using vectors instead of the square 
matrix iJj, with possible advantages in storage and number of operations. For details 
see Abaffy and Spcdicato (1989) and on how these formulations can be used to imbed 
in an ABS approach even iterative methods (including the Kaczmarz, Gauss-Seidel, 
ORTHOMIN, ORTHODIR etc. methods), see Spedicato and Li (2000). 

The choice of the parameters ifi, Vi, ,Zi, Wi determines particular methods. 
The basic ABS class is obtained by taking Vi = Ci, as the i-th unit vector in K^. 

We recall some properties of the scaled ABS class, assuming that A has full rank. 

• Define Vi — {vi, . . . ,Vi) and Wi — {wi, . . . , Wi). Then 

Hi+,A'^Vi = 0, Hl,,Wi = (8) 



• The vectors HiA^Vi, Hjwi are zero if and only if a^, Wi are respectively linearly 
dependent from oi, . . . , aj_i, Wi, . . . , Wi^\. 
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• Define Pi — (pi, ■ ■ ■ ,Pi) and Ai — {ai, . . . ,ai). Then the imphcit factorization 
V^AfPi = Li holds, where Li is nonsingular lower triangular. Hence, if m = n, 
one obtains a semiexplicit factorization of the inverse, with P — Pn, V — 

Kl; L = Ln 

A-^ = PL-^V^. (9) 

For several choices of V the matrix L is diagonal, hence formula (8) gives a 
fully explicit factorization of the inverse as a byproduct of the ABS solution 
of a linear system. 

• The general solution of system (1) can be written as follows, with g e i?" 
arbitrary 

X = Xm+l + H^+iQ (10) 

• The Abaffian can be written in terms of the parameter matrices as 

Hi+, ^H,- H,A'^Vi{W^H,A^Vi)-'W^H^. (11) 

Letting V = Vm, W = Wm, one can show that the parameter matrices Hi, V, W are 
admissible (i.e. condition (7) is satisfied) iff the matrix Q = V^AHfW is strongly 
nonsingular (i.e. it is LU factorizable) . This condition can be satisfied by exchanges 
of the columns of F or 1^. If Q is strongly nonsingular and we take, as is done in all 
algorithms so far considered, Zi = Wi, then condition (4) is also satisfied. Analysis 
of the conditions under which Q is not strongly nonsingular leads, when dealing 
with Krylov space methods in their ABS formulation, to a characterization of the 
topology of the starting points that can produce a breakdown (either a division by 
zero or a vanishing search direction) and to several ways of curing it, including those 
considered in the literature. 

Three subclasses of the scaled ABS class and particular algorithms are now 
recalled. 

(a) The conjugate direction subclass. This class is obtained by setting t',; = pi. It 

is well defined under the condition (sufficient but not necessary) that A is 
symmetric and positive definite. It contains the ABS versions of the Choleski, 
the Hestenes-Stiefel and the Lanczos algorithms. This class generates all pos- 
sible algorithms whose search directions are 74-conjugate. If Xi — 0, the vec- 
tor minimizes the energy (A-weighted Euclidean) norm of the error over 
Span{pi, ■ ■ ■ ,Pi) and the solution is approached monotonically from below in 
the energy norm. 

(b) The orthogonally scaled subclass. This class is obtained by setting Vi = Api. 

It is well defined if A has full column rank and remains well defined even if 
m is greater than n. It contains the ABS formulation of the QR algorithm 
(the so called implicit QR algorithm), the GMRES and the conjugate residual 



5 



algorithms. The scahng vectors are orthogonal and the search vectors are A^A- 
conjugate. If xi = 0, the vector Xj+i minimizes the Euclidean norm of the 
residual over Span{pi, . . . ,pi) and the solution is monotonically approached 
from below in the residual norm. It can be shown that the methods in this 
class can be applied to overdetermined systems of m > n equations, where in 
n steps they obtain the solution in the least squares sense. 

(c) The optimally stable subclass. This class is obtained by setting Vi = {A~^)-'^pi, 
the inverse disappearing in the actual recursions. The search vectors in this 
class are orthogonal. If Xi — 0, then the vector Xi+i is the vector of least 
Euclidean norm over Span{pi, ■ ■ ■ ,Pi) and the solution is approached mono- 
tonically from below in the Euclidean norm. The methods of Gram-Schmidt 
and of Craig belong to this subclass. The methods in this class have minimum 
error growth in the approximation to the solution according to a criterion by 
Broyden. 



3 Special algorithms in the basic ABS class 

In this section we consider in more detail three specially important algorithms, that 
belong to the basic ABS class (i.e. V = I). 

(A) The Huang algorithm is obtained by the choices Hi = I, Zi = Wi = ai. A 
mathematically equivalent, but numerically more stable, formulation is the so 
called modified Huang algorithm {pi = Hi^Hitti) and i^j+i = H^ — PipJ /pjpi)- 
Huang algorithm generates search vectors that are orthogonal and identical 
with those obtained by the Gram-Schmidt procedure applied to the rows of 
A. If xi = 0, then x^+i is the solution with least Euclidean norm of the first i 
equations. The solution x^ with least Euclidean norm of the whole system is 
approached monotonically and from below by the sequence Xi. For arbitrary 
starting point the Huang algorithm generates that solution which is closest in 
Euclidean norm to the initial point. 

(B) The implicit LU algorithm is given by the choices Hi = I, Zi = Wi = e^. It 
is well defined iff A is regular (i.e. all principal submatrices are nonsingu- 
lar). Otherwise column pivoting has to be performed (or, if m = n, equation 
pivoting) . The Abaffian has the following structure, with Ki e 



H 











(12) 



implying that the matrix Pj is unit upper triangular, so that the implicit fac- 
torization A — LP~^ is of the LU type, with units on the diagonal. The 



6 



algorithm requires ior m — n, multiplications plus lower order terms, 

the same cost of classical LU factorization or Gaussian elimination. Storage 
requirement for Ki requires at most positions, i.e. half the storage needed 
by Gaussian elimination and a fourth that needed by the LU factorization algo- 
rithm (assuming that A is not overwritten). Hence the implicit LU algorithm 
has same arithmetic cost but uses less memory than the most efficient classical 
methods. 

(C) The implicit LX algorithm, see Spedicato, Xia and Zhang (1997), is defined 
by the choices Hi — I, Zi — Wi — Cfc., where ki is an integer, 1 < ki < n, 
such that e^Hitti 7^ 0. By a general property of the ABS class for A with full 
rank there is at least one index ki such that e^Hiai 7^ 0. For stability reasons 
we select ki such that | e^Mitti \ is maximized. This algorithm has the same 
overhead and memory requirement as the implicit LU algorithm, but does not 
require pivoting. Its computational performance is also superior and generally 
better than the performance of the classical LU factorization algorithm with 
row pivoting, as available for instance in LAPACK or MATLAB, sec Mirnia 
(1996). Therefore this algorithm can be considered as the most efficient general 
purpose linear solver not of the Strassen type. 

4 Solution of linear Diophantine equations 

One of the main results in the ABS field has been the derivation of ABS methods 
for linear Diophantine equations. The ABS algorithm determines if the Diophan- 
tine system has an integer solution, computes a particular solution and provides a 
representation of all integer solutions. It is a generalization of a method proposed 
by Egervary (1955) for the particular case of a homogeneous system. 

Let Z be the set of all integers and consider the Diophantine hnear system of 
equations 

Ax = b, xeZ"", Ae Z™^", beZ"", m<n. (13) 

While thousands of papers have been written concerning nonlinear, usually poly- 
nomial, Diophantine equations in few variables, the general linear system has at- 
tracted much less attention. The single linear equation in n variables was first solved 
by Bertrand and Betti (1850). Egervary was probably the first author deahng with 
a system (albeit only the homogeneous one). Several methods for the nonhomoge- 
neous system have recently been proposed based mainly on reduction to canonical 
forms. 

We recall some results from number theory. Let a and h be integers. If there is 
integer 7 so that b — then we say that a divides b and write a\b, otherwise we 

write aJI b. If ai, . . . , a„ arc integers, not all being zero, then the greatest common 
divisor (gcd) of these numbers is the greatest positive integer 6 which divides all 
Oj, i = 1, . . . ,n and we write 6 = gcd{ai, . . . , a„). We note that 5 > 1 and that 
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5 can be written as an integer linear combination of the a^, i.e. 5 — z^a for some 
z e i?". One can show that 5 is the least positive integer for which the equation 
aiXi + ... + QnXn — S has an integer solution. Now S plays a main role in the 
following 

Fundamental Theorem of the Linectr Diophantine Equation 

Let Oi, . . . , a„ and b be integer numbers. Then the Diophantine linear equation aiXi + 
. . . + a„a;„ = b has integer solutions if and only if gcd{ai, . . . , a„)| b. In such a case 
if n > 1 then there is an infinite number of integer solutions. 

In order to find the general integer solution of the Diophantine equation aiXi + . . . + 
dnXn — b, the main step is to solve aiXi + . . . + a„a;„ = S, where S = gcd{ai, . . . , a„), 
for a special integer solution. There exist several algorithms for this problem. The 
basic step is the computation of 6 and z, often done using the algorithm of Rosser 
(1941), which avoids a too rapid growth of the intermediate integers, and which 
terminates in polynomial time, as shown by Schrijver (1986). The scaled ABS algo- 
rithm can be applied to Diophantine equations via a special choice of its parameters, 
originating from the following considerations and Theorems. 

Suppose Xi is an integer vector. Since Xj+i = Xi — aiPi, then Xj+i is integer 
if ai and pi are integers. If vfApi\{vfri), then is an integer. If Hi and Zi are 
respectively an integer matrix and an integer vector, then pi — Hfzi is also an 
integer vector. Assume Hi is an integer matrix. Prom (6), if vj AHfwi divides all 
the components of HiA^Vi, then Hi^i is an integer matrix. 

Conditions for the existence of an integer solution and determination of all integer 
solutions of the Diophantine system are given in the following theorems, generalizing 
the Fundamental Theorem, see Esmaeih, Mahdavi-Amiri and Spedicato (2001a). 

Theorem 1 Let A be full rank and suppose that the Diophantine system (13) 
is integer solvable. Consider the Abaffians generated by the scaled ABS algorithm 
with the parameter choices: Hi is unimodular (i.e. both Hi and Hi^ are integer 
matrices); for i = 1, . . . ,m, Wi is such that wfHiA^Vi — Si, Si — gcd{HiA^Vi). 
Then the following properties are true: 

(a) the Abaffians generated by the algorithm are well-defined and are integer matri- 

ces 

(b) ify is a special integer solution of the first i equations, then any integer solution 

X of such equations can be written as x — y-\- H^^^q for some integer vector q. 

Theorem 2 Let A be full rank and consider the sequence of matrices Hi gen- 
erated by the scaled ABS algorithm with parameter choices as in Theorem 1. Let Xi 
in the scaled ABS algorithm be an arbitrary integer vector and let Zi be such that 
zfHiA^Vi = gcd{HiA^Vi). Then system (13) has integer solutions iff gcd{HiA^Vi) 
divides vjri for i — 1, . . . ,m. 
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Prom the above theorems we obtain the following scaled ABS algorithm for Dio- 
phantine equations. 

The ABS Algorithm for Diophantine Linear Equations 

(1) Choose Xi e Z"-, arbitrary, Hi e Z"-^"-, arbitrary unimodular. Let i = 1. 

(2) Compute Tj = vfri and Sj = HiA^Vi. 

(3) // (sj = and Tj = 0) then let Xj+i = Xi, i^^+i = H^, rj+i = rj and go to step 
(5) (the ith equation is redundant). If {si = and Tj ^ 0) then Stop (the ith 
equation and hence the system is incompatible). 

(4) {sj 7^ 0} Compute 5i — gcd{si) and Pi — Hfzi, where Zi e is an arbitrary 
integer vector satisfying zjsi — 5i. If SiJ(Ti then Stop (the system is integerly 
inconsistent), else Compute ccj = Ti/Si, let Xj+i = Xi — aipi and update Hi by 
Hi^i = Hi — ^^ ' h-a^v^' where Wi e is an arbitrary integer vector satisfying 
wfsi = 5i. 

(5) Ifi — m then Stop (x^+i is a solution) else let i — i + 1 and ^'o to step (2). 

It follows from Theorem 1 that if there exists a solution for the system (13), then 
X — Xm+i + H^_^iq, with arbitrary q e Z", provides all solutions of (13). 

Egervary's algorithm for homogeneous Diophantine systems corresponds to the 
choices Hi = I, Xi = and Wi = Zi, for all i. Egervary claimed, without proof, 
that any set of n — m linearly independent rows of i?rre+i form an integer basis 
for the general solution of the system. Esmaeih et al. (2001a) have shown by a 
counterexample that Egervary's claim is not true in general; we have also provided 
an analysis of conditions under which m rows in Hm+i can be eliminated. 

One can show that by special choices of the scaling parameters and by an inessen- 
tial modification of the update of the Abaffian it is not necessary to solve in general 
2m single n-dimensional Diophantine equations (to determine Zi, Wi), but just one 
single Diophantine equation, at the last step. 

The ABS algorithm for linear Diophantine systems has been extended also to 
systems of m linear inequalities (m < n), where it provides an almost explicit repre- 
sentation of all solutions, see Esmaeili, Mahdavi-Amiri and Spedicato (2001b). The 
used technique allows, in the case of m < n real inequalities to find a fully explicit 
representation of all solutions, see Esmaeili, Mahdavi-Amiri and Spedicato (2000), 
thereby obtaining a generalization of formula (10). Finally, the ABS approach pro- 
vides a new way of computing, given an integer matrix, its integer null space and, 
under some conditions, an integer basis of it, see Esmaeili, Mahdavi-Amiri and 
Spedicato (2001c). 
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5 Reformulation of the simplex method via the 
implicit LX method 

The implicit LX algorithm has a natural application in reformulating the simplex 
method for the LP problem in standard form: 

min (Fx 

subject to 

Ax = h, x>Q. 

The suitability of the implicit LX algorithm derives from the fact that the al- 
gorithm, when started with xi = 0, generates basic type vectors Xj+i, which are 
vertices of the polytope defined by the constraints of the LP problem if the compo- 
nents not identically zero are nonnegative. 

The basic step of the simplex method is moving from one vertex to another 
according to certain rules and reducing at each step the value of Fx. The direction 
of movement is obtained by solving a linear system, whose coefficient matrix Ab, 
the basic matrix, is defined by m linearly independent columns of A, called the basic 
columns. This system is usually solved via the LU factorization method or also 
via the QR factorization method. The new vertex is associated with a new basic 
matrix A'^, obtained by substituting one of the columns ol Ab with a column of 
the matrix Aat that comprises the columns of A not belonging to Ab. One has 
then to solve a new system where just one column has been changed. If the LU 
factorization method is used, then the most efficient way to recompute the modified 
factors is the steepest edge method of Forrest and Goldfarb (1992), that requires 

multiplications. This implementation of the simplex method has the following 
storage requirement: m? for the LU factors and mn for the matrix A, that has to 
be kept to provide the columns needed for the exchanges. 

The reformulation of the simplex method via the implicit LX method, see Zhang 
and Xia (1995) and Spedicato, Xia and Zhang (1997), exploits the fact that in the 
implicit LX method the exchange of the j-th column m Ab with the k-ih. column 
in Ajv corresponds to exchanging previously chosen parameters Zj — Wj = Cj^ 
with new parameters Zk — — e^g. in terms of the implied modification of the 
Abaffian this operation is a special case of a general rank-one modification of the 
parameter matrix W = {wk^, . . . ,Wk^). The modified Abaffian can be efficiently 
evaluated using a general formula of Zhang (1995), without explicit use of the k-th 
column oi An. Moreover all information needed to build the search direction (the 
polytope edge) and to implement the (implicit) column exchange is contained in 
the Abaffian matrix. Thus there is no need to keep the matrix A. Hence storage 
requirement is only that needed in the construction and keeping of the matrix //m+i; 
i.e. respectively and n(n — m). For m close to n, such storage is about 8 times 
less than in the Forrest and Goldfarb implementation of the LU method. For small m 
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storage is higher, but there is an alternative formulation of the implicit LX method 
that has a similar storage, see Spedicato and Xia (1999). 

We now give the main formulas for the simplex method in the classical and in 
the ABS formulation. The column in An that substitutes a column in ^4^ is usually 
taken as the column with least relative cost. In the ABS approach this corresponds 
to minimize with respect to i e N^n the scalar rji — c^H^Ci. Let N* be the index 
so chosen. The column in to be exchanged is often chosen with the criterion 
of the least edge displacement that keeps the basic variables nonnegative. Defining 
LUi = x^Ci/eJ H'^cn* with x the current vertex, the above criterion is equivalent to 
minimize uji with respect to the set of indices i G Bm such that 

ejH^eN* > (14) 

Notice that H'^cn* 7^ and that an index i such that (14) be satisfied always exists, 
unless X is a solution. 

The update of the Abaffian after the exchange of the obtained unit vectors is 
given by the following simple rank-one correction 

H'^H- (Hcb* - eB*)e^.Hle^.HeB^ (15) 

The displacement vector classically obtained by solving system Asd = —Ae^*, is 
obtained at no cost hy d — H^_^^eN*- The relative cost vector, classically given by 
r = c—A^Aq^cb, with cb comprising the components of c with indices corresponding 
to those of the basic columns, is given by formula r = Hj^+ic. 

It is easily seen that update (15) requires no more than m{n—m) multiplications. 
Cost is highest for m = n/2 and gets smaller as m gets smaller or closer to n. In the 
method of Forrest and Goldfarb multiplications are needed to update the LU 
factors and again to solve the system. Hence the ABS approach via formula (15) 
is faster for m > n/3. For m < n/3 similar costs are obtained using the alternative 
formulation of the implicit LX algorithm described in Spedicato and Xia (1999). For 
m very close to n the advantage of the ABS formulation is essentially of one order 
and such that no need is seen to develop formulas for the sparse case. 

There is an ABS generalization of the simplex method based upon a modification 
of the Huang algorithm which is started by a certain singular matrix, see Zhang 
(1997). In such generalization the solution is approached via a sequence of points 
lying on the faces of the polytope. If one of such points happens to be a vertex, then 
all successive iterates are vertices and the standard simplex method is reobtained. 

6 ABS unification of feasible direction methods 
for linearly constrained minimization 

ABS methods allow a unification of feasible direction methods for linearly con- 
strained minimization, of which the LP problem is a special case. Let us first 
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consider the problem with only equality constraints 

min f{x), x e i?" 

subject to 

Ax^b, Ae -R""'", m<n, rank{A) = m. 

Let Xi be a feasible initial point. If we consider an iteration of the form Xj+i = 
Xi — aidi, the sequence Xi consists of feasible points iff 

Adi = (16) 

The general solution of (16) can be written using the ABS formula (10) 

di = Hl^.q (17) 

In (17) matrix Hm+i depends on parameters Hi, W and V and q G i?" can also be 
seen as a parameter. Hence the general iteration generating feasible points is 

Xi+i ^Xi- aiH'^_^-^q. (18) 

The search vector is a descent direction if d^Vf{x) = q^ Hm+i^ f{x) > 0. This 
condition can always be satisfied by a choice of q unless iJ^+i V/(a:) = 0, implying, 
from the structure of NuU{H^^i), that V/(x) = A^\ for some A, hence that Xj+i 
is a KT point with A the vector of Lagrange multipliers. If Xj+i is not a KT point, 
we can generate descent directions by taking 

q = QHm+iVf{x) (19) 

with Q symmetric positive definite. We obtain therefore a large class of methods 
with four parameter matrices {Hi,W,V,Q). 

Some well-known methods in the hterature correspond to taking W — I in (19) 
and building the Abaffian as follows. 

(a) The reduced gradient method of Wolfe. Hm+i is built via the implicit LU 

method. 

(b) The projection method of Rosen, ifm+i is built via the Huang method. 

(c) The Goldfarb and Idnani method. Hm+i is built via a modification of Huang 

method where Hi is a symmetric positive definite approximation of the inverse 
Hessian of f{x). 

To deal with linear inequality constraints there are two approaches in literature. 
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• The active set method. Here the set of equahty constraints is augmented with 

some inequaUty constraints, whose selection varies in the course of the process 
till the final set of active constraints is determined. Adding or cancelling a 
single constraint corresponds to a rank-one correction to the matrix defining 
the active set. The corresponding change in the Abaffian can be performed in 
order two operations, see Zhang (1995). 

• The standard form approach. Here one uses slack variables to put the problem 
in the following equivalent standard form 

min f{x), 

with the constraints 

Ax ^b, X >0. 

If Xi satisfies the above constraints, then a sequence of feasible points is generated 
by the iteration 

Xi+i ^Xi- aiPiHm+iVf{x) (20) 

where CKj may be chosen by a line search along vector Hm+i'V f{x), while /3i > is 
chosen to avoid violation of the nonnegativity constraints. 

If f{x) is nonlinear, then Hm+i can usually be determined once for all at the first 
iteration, since generally V f{x) changes from point to point, allowing the determi- 
nation of a new search direction. However if f[x) = c^x is linear, in which case we 
obtain the LP problem, then to get a new search direction we have to change Hm+i- 
We already observed in section 4 that the simplex method in the ABS formulation 
is obtained by constructing the matrix Hm+i via the implicit LX method and at 
each step modifying one of the unit vectors used to build the Abafiian. One can 
show that the Karmarkar method corresponds to Abaffians built via a variation of 
the Huang algorithm, where the initial matrix is Hi — Diag{xi) and is changed 
at every iteration (whether the update of the Abaffian can be performed in order 
two operations is still an open question). We expect that better methods may be 
obtained by exploiting all parameters available in the ABS class. 

7 ABSPACK and its numerical performance 

The ABSPACK project (based upon a collaboration between the University of Berg- 
amo, the Dalian University of Technology and the Czech Academy of Sciences) aims 
at producing a mathematical package for solving linear and nonlinear systems and 
optimization problems using the ABS algorithms. The project will take several years 
for completion, in view of the substantial work needed to test the alternative ways 
ABS methods can be implemented (via different linear algebra formulations of the 
process, different possibilities of reprojections, different possible block formulations 
etc.) and of the necessity of comparing the produced software with the estabhshed 
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packages in the market (e.g. MATLAB, LINPACK, LAPACK, UFO ...). It is ex- 
pected that the software will be documented in a forthcoming monograph and wilU 

be made available to general users. 

Presently FORTRAN 77 implementations have been made of several versions of 
the following ABS algorithms for solving linear systems: 

1. The Huang and the modified Huang algorithms in two different linear algebra 
versions of the process, for solving determined, underdetermined and overde- 
termined systems, for a solution of least Euclidean norm 

2. The implicit LU and implicit LX algorithm for determined, underdetermined 
and overdetermined linear systems, for a solution of basic type 

3. The implicit QR algorithm for determined, underdetermined and overdeter- 
mined linear systems, for a solution of basic type 

4. The above algorithms for some structured problems, namely KT equations 
and banded matrices. 

For a full presentation of the above methods and their comparison with NAG, LA- 
PACK, LINPACK and UFO codes see Bodon, Spcdicato and Luksan(2000a,b,c). 
Some results are presented at the end of this section. There the columns refer re- 
spectively to: the problem, the dimension, the algorithm, the relative solution error 
(in Euclidean norm), the relative residual error in Euclidean norm (i.e. ratio of resid- 
ual norm over norm of right hand side), the computed rank and the time in seconds. 
Computations have been performed in double precision on a Digital Alpha worksta- 
tion with machine zero about 10^^^. All test problems have been generated with 
integer entries or powers of two such that all entries are exactly represented in the 
machine and the right hand side can be computed exactly, so that the given solution 
is an exact solution of the problem as it is represented in the machine. Compari- 
son is given with some LAPACK codes, including those based upon singular value 
decomposition (svd) and rank revealing QR factorization {gqr). 
Analysis of all obtained results indicates: 

1. Modified Huang is generally the most accurate ABS algorithm and compares 
in accuracy with the best LAPACK solvers based upon singular value fac- 
torization and rank revealing QR factorization; also the estimated ranks are 
usually the same. 

2. On problems where the numerical estimated rank is much less than the dimen- 
sion, one of the versions of modified Huang is much faster than the LAPACK 
codes using SVD or rank revealing QR factorization, by even a factor more 
than 100. This is due to the fact that once an equation is recognized as de- 
pendent it does not contribute to the general overhead in ABS algorithms. 
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3. Modified Huang is more accurate than other ABS methods and faster than 
classical methods on KT equations. 

It should be noted that the performance of the considered ABS algorithms in term 
of time could be improved by developing block versions. ABS methods are moreover 
expected to be faster than many classical methods on vector and parallel computers, 
see Bodon (1993). 
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Numerical Results 



RESULTS ON DETERMINED LINEAR SYSTEMS 
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RESULTS ON OVERDETERMINED SYSTEMS 



Condition number: 0.16D+21 
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RESULTS ON UNDERDETERMINED LINEAR SYSTEMS 



Condition number: 0.29D+18 
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RESULTS ON KT SYSTEMS 



Condition number: 0.26D+21 
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